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A temporal study of energy transfer across length scales is performed in 3D numer¬ 
ical simulations of homogeneous shear flow and isotropic turbulence. The average 
time taken by perturbations in the energy flux to travel between scales is measured 
and shown to be additive. Our data suggests that the propagation of disturbances 
in the energy flux is independent of the forcing and that it defines a ‘velocity’ that 
determines the energy flux itself. These results support that the cascade is, on 
average, a scale-local process where energy is continuously transmitted from one 
scale to the next in order of decreasing size. 


The difficulty in understanding a multiscale problem such as turbulence has been fre¬ 
quently tackled by assuming a priori some type of simplified phenomenology. Take, for 
instance, Richardson’s cartoon based on concepts such as eddies whose energy cascades by 
the successive breakup of larger eddies into smaller ones, until it is dissipated by viscosity. 1 
It was later used as the basis for the more quantitative work of Kolmogorov. We now know 
that in 3D turbulence the energy does cascade towards the smallest scales, ’’ at least on 
average. But several models have been discussed in the literature which are consistent with 
this average trend yet differ in their detailed mechanism. For example, the energy could 
jump directly from large eddies to much smaller ones, 1 ’ ’ contradicting the scale locality 
assumed by Richardson, or include frequent excursions from smaller to larger scales - a 
process coined backscatter - questioning a unique directionality. ’ Such alternative roads 
leading to the same average behavior urge an improved understanding of the cascade dy¬ 
namics. In this letter, we report on the time taken by disturbances in the the energy flux 
to travel in scale space - an essential ingredient in unsteady phenomenological models. 
To overcome the difficulty in generating a wide dynamic range in space but especially in 
time, our direct numerical simulations are purposely run for very long times to perform a 
temporal cross-correlation analysis between energy fluxes at various length scales. 

Our approach is based on the large- and small-scale decomposition of the instantaneous 
velocity field according to Ui(xi,t) = Ui(xi,t) + u'^Xi, t),' where Ui is the i-th component of 
the spatially low-pass filtered velocity. In an incompressible flow with kinematic viscosity 
v, the kinetic energy of the large-scale field evolves as 

/ <9 ~<9\l_ d ~ \ ~ ~ ~ ~ 

T Uj — — J —Uilli — d - IliTij QlSUiSij J 2lS Sij Sij E T Uifi, (1) 

where Sij = ( dui/dxj + dv,j/dxi) /2 is the strain-rate tensor of the large-scale velocities, 
E = — TijSij , Tij = ufuj — UiUj is the subgrid-scale stress tensor and fi is the forcing term. 
Eq. (1) has been studied extensively in the context of large-eddy simulations (LES), where 
capturing the behavior of E is at the cornerstone of most modeling strategies.' 1 In the mean, 
E is positive and acts as a net energy removal from the large scales by the small ones. 1,1 
We choose it as our real-space marker of cross-scale energy transfer, and study it in two 
different flows: homogeneous shear turbulence (HST) and homogeneous isotropic turbulence 
(HIT). The details of the simulations are given in Table I. The low-pass filtered velocity 
in HIT was obtained by multiplying the Fourier modes Ui{k,t) by an isotropic Gaussian 
kernel G(k) = exp[— (rfc) 2 /24], 11 ’ 1 where k is the wavevector and r the filter width. In 


djose@torroj a. dmt .upm.es 



2 


TABLE I. Parameters of the simulations. Rex is the Reynolds number based on the Taylor- 
microscale. Ni and Li are the number of real Fourier modes and the domain size in directions 
i = x,y,z. Length scales are p = (i' 3 /e) 1,/4 and L 0 = K 3 ^ 2 /e. Times are normalized by T 0 = K/e. 
Tsimu is the simulation time, and A tt.ot is the average delay between (K) and (e). Tkk is an 
autocorrelation time for (K), defined in the text. K = mUi/2. 


Case 

Rex 

N x xN y xN z 

( L x x Ly x L z )/rj 

Lo/v 

Tsimu/To 

A ttot/To 

Tkk /Aftot 

HST 

107 

768 x 512 x 255 

1117 x 745 x 372 

267 

213 

0.49 

9.04 

HIT1 

146 

256 3 

506 3 

425 

165 

0.53 

2.32 

HIT2 

236 

512 3 

1011 3 

876 

16 

0.40 

2.36 

HIT3 

384 

1024 3 

2022 3 

1813 

4 

0.42 

2.14 


HST this method could only be applied in the streamwise and spanwise directions, along 
which the flow was Fourier-discretized. Since seventh-order compact finite differences were 
used in the vertical (sheared) direction, filtering along it was implemented by convolving 
the velocity with the real-space transform of G (k ): G (x) = P exp[—6x 2 /r 2 ], where P is 
a constant chosen to meet the normalization condition. For comparison, a sharp spectral 
filter was used on some HIT fields with G(k) = 0 when |fc| > tv/ r, G(k) = 1 otherwise. 

Other energy transfer markers exist which could have been suited to the F1IT simulation. 
An obvious candidate comes from the spectral energy equation in isotropic turbulence 

^^-=F(k,t)-2uk 2 E(k,t) + E(k,t), (2) 

where E ( k,t ) is the 3D instantaneous energy spectrum, k = |fc| and S is the forcing. The 
spectral energy flux II(fc) = F(k)dk is often invoked in energy cascade studies. For a 
single flow field of isotropic turbulence, n(fc) and (£(r)} are equal when a sharp spectral 
filter is used to compute £ with a cut-off wavenumber k = tv/t. (Hereafter, (9) is the 
time-dependent spatial average of 6 over the computational domain, while 9 is the temporal 
mean of (9).) The difference between the two energy transfer markers thus amounts to 
a choice of filter type. Since n is of limited use in flows other than HIT and given the 
relevance of £ in LES, we favored £ for comparison between the two flows. 

The probability density function p(£) can be seen on Fig. 1(a). Its positive skewness 
indicates that strong events are more likely when £ acts as a sink in Eq. (1) than as a 
source. The tails of p become narrower with increasing filter width, in agreement with 
Ref. 11. Whereas the volume ratio of forward cascade to backscatter is known to favor 
the former over the latter for Gaussian filters/ 5 ’ 1 the ratio of the total forward energy flux 
>1 >f = / 0 °° £p (£) d£ to that of backscatter &b = f £p (£) d£ is not documented. Their 
ratio is given on Table II, which reveals an even stronger preponderance of the forward 
cascade than what could be inferred from the volume ratios alone. To emphasize this point, 
Fig. 1(c) displays the value of the integrand defining and The sharp spectral filter 
leads to a much more symmetric picture of the cascade, as inferred by comparing Fig. 1(c) 
and (d). Fig. 1(b), with a sharp filter, shows that p(£) is less skewed than on Fig. 1(a), 
with a Gaussian filter, yet the wider tails for decreasing r occur with both filters. 


TABLE II. Ratio $f/$s of forward to reverse energy flux. Numbers in parentheses are the 
corresponding volume ratios, J 0 °° p (£) d£/ f_ p (£) d£ . 


r/v 

8 

16 

31 

62 

125 

250 

HST(Gauss) 

N/A 

16(4) 

24(5) 

62(10) 

139(19) 

N/A 

HIT2 (Gauss) 

15(3) 

15(4) 

21(5) 

39(8) 

59(10) 

37(9) 

HIT2 (sharp) 

1.4(1.2) 

1.7(1.3) 

2.2(1.6) 

2.9(1.9) 

3.7(2.3) 

5.0(2.7) 
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FIG. 1. Top: probability density functions (PDFs) of E in HIT2, normalised by the standard 
deviation E\ (a) Gaussian filter, (b) Sharp filter. Bottom: weighted PDF, whose integral defines 
and 4 >b. (c) Gaussian filter, (d) Sharp filter. 


We now move towards the dynamics of (E), based on its time series which we computed 
for all flows with a Gaussian filter - the sharp filter was used only on a few HIT fields 
widely spaced in time. They are shown on Fig. 2(a) at two filter widths, together with 
the kinetic energy ( K) = ( mui ) /2 and the viscous dissipation (e). It is clear that the 
signals are correlated, but that there is a delay separating them. ( K) and (e) behave 
as the earliest and latest signals, while the delay of (E) with respect to (K) increases 
with decreasing r. To visualize this effect more clearly, Fig. 2(c) displays the temporal 
evolution of (E) as a color map where the abscissae are time. Color bands corresponding 
to (E) have been ordered vertically with r decreasing logarithmically downwards, and (e) 
added at the bottom. The propagation across r and t of disturbances in (E) is evident. 
To quantify this process, we compute the temporal cross-correlation of all these signals 
with each other, as well as their temporal autocorrelation. A few such correlations are 
illustrated on Fig. 2(b), where the peak appears at the average delay between the two 
chosen signals. We start by looking at Ai tot , the average time taken by a change in (. K) 
to propagate and appear as a change in (e). It is compiled in Table I for all our flows, 
which shows that At tot is approximately half the integral dissipation time T 0 = K/e. 
Similar data are scarce in the literature, so that a comparison is not straightforward. In a 
computational study of homogeneous shear flow, the delay between the time histories of 
(I\) and (e) was estimated to be of the order of A t t0 t/T 0 « 1.5 at Re\ « 50. Since this 
disagrees with our A t t0t /T 0 = 0.49 at Re\ ~ 107, we repeated the simulation in Ref. 14, 
and found At to t/T 0 = 0.52. The discrepancy can probably be attributed to the different 
estimation methods. Cross-correlations were not computed in Ref. 14, where A ttot was 
not the main focus of the study. A value of Att. 0 t/T 0 = 0.44 can be extracted from the 
data in Ref. 15 where a DNS of HIT was ran at Re\ = 122, which agrees well with our 
findings. In Ref. 16, a A ttot was defined as the delay between (K) 3 ^ 2 /Li n t and (e), where 
Lint = (37t/4(A')) J k~ 1 E(k) dk. Using temporal cross-correlations in HIT at Re\ = 219 
- comparable to our HIT2, they found Af t0t /T 0 ps 0.21. We found At t0t /T 0 ps 0.28 using 
our HIT2 and the same definition of A ttot as in Ref. 16, confirming that the lag between 
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FIG. 2. (a) Temporal evolution of spatially-averaged quantities, centered and normalized by their 
standard deviation; HIT2. (b) Cross-correlation curves between time series of (E) at various filter 
widths and between (K) and (e); HIT2. (c) Time-scale diagram of (E), with r decreasing from top 
to bottom and (e) added as the bottom band; HIT3 with r/r) values from Table II. The dash-dotted 
line corresponds to e 1 ,/ 3 A t = (250 r/) 2 ^ 3 — r 2 / 3 - see Eq. ( 6 ). 


(K ) 3 / 2 /L int and (e) is shorter than that between ( K) and (e). 

A comment on the origin of the time dependence of ( K ), (e) and (£) is in order. The 
temporal oscillations of (K) and ( e ) in HST are known to be physically caused and related 
to bursting . In the HIT simulations, the turbulence is sustained by a deterministic force 


fi (M) 


£Ui (fc, t) / [2 Ef(t)\ , if 0 < k < kf, 
0, otherwise, 


( 3 ) 


where e is the target mean dissipation, Ef(t ) = f^ f E{k,t) d k and kf = Att/L x . This 
commonly used scheme is mildly unstable, because of the delay between fi and (e). It 
generates time oscillations of the energy while maintaining a constant resolution of k max r] = 
1.5 in the mean. With these two flows at hand, we can assess the dependence of Ai tot on 
the large-scale forcing. We define the characteristic time scale Tkk of the kinetic energy 
as the width of the temporal autocorrelation of ( K) at half its peak height. The ratio 
Tkk/A ttot is between 2 and 2.5 for all our HIT simulations, but about 9 for the HST - 
see Table I. Yet changes in ( K) appear as changes in (e) within half a large-eddy turnover 
time in the two differently forced flows, suggesting that A tt 0 t/T 0 is a common feature of the 
energy cascade when the large scales fluctuate with periods in our range of Tkk/T 0 . The 
dependence of our measured A t to t/T 0 on the large-scale period could be studied further by 
extending this range of Tkk/T a with the addition of a modulating frequency in the forcing, 
as done Ref. 18, or with a stochastic forcing. 


TABLE III. Symbol legend for Figs. 3 and 5. r a /r] = 10 ay/E/n so that r 1 m 8 r/, r 2 ~ 1677 , etc. 



r (e> 

a = 1 

a = 2 

a = 4 

a = 8 

a = 16 

HST 

+ 

N/A 

X 

* 

O 

N/A 

HIT1 

□ 

A 

0 

V 

N/A 

N/A 

HIT2 

■ 

▲ 

• 

T 

► 

N/A 

HIT3 

■ 

A 

• 

T 

► 

< 
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FIG. 3. (a) Delay between (S) and (e), measured either as the one-step delay, or as the sum of 
two intermediate steps involving r a . (b) Ratio between two- and one-step delays. See Table III for 
symbols of r a and flow. 


We next test the additivity of the delay times. We want to see if the time taken by a 
disturbance in the energy flux in going from scale size r 3 to r\ is equal to the sum of the 
delays in going from r 3 to r 2 and from r 2 to rq, where > r 2 > ri- Fig. 3(a) displays the 
average time needed for disturbances in the energy fluxes at a given scale r to travel down 
to (e). This time is computed for all available combinations of two intermediate jumps 
starting at r and ending in (e). The agreement between the different jump combinations 
within the same flow and across different flows is satisfactory. In order to highlight any 
residual discrepancy, we display the value of the ratios between one- and two-step cascading 
times on Fig. 3(b). A ratio close to unity implies additivity of the delays, which is confirmed 
for all the starting scales r and flows examined. Note the narrow range of the vertical axis. 
Note also that Fig. 3(a) hints at two different regimes for r above and below approximately 
30fp This is consistent with the results of Ref. 20, who showed that viscous eddies below 
r/p ~ 30 are enslaved to larger ones above that scale. In essence, r/p = 30 is the lower 
limit of the inertial cascade. 

In an influential paper, Lumley discusses two cascade models, 1 and proposes a forcing 
experiment similar to the present one to distinguish between them. He starts by considering 
a hierarchy of discrete eddies of decreasing size. In the first model, each eddy transfers its 
energy to the one immediately below. The transfer occurs at a rate determined by the 
corresponding scale-dependent eddy turnover time, [k 3 E(k)]~ 1 ^ 2 ~ k ~ 2 / 3 ~ r 2 / 3 , so that 
the propagation of energy from the large towards the small scales develops into a front-like 
diffusion through scale-space with a finite scale-dependent velocity. In the second model, 
most of the energy is still transferred to the immediately smaller eddy below, but a fraction 
is passed to other eddies further along the cascade. Hence the smallest eddies receive a small 
amount of energy almost immediately after it is injected into the system, and increasing 
amounts as time goes on. The difference between the two models is that all the energy 
in the first one has to pass through each eddy size, resulting in additive cascade times, 
while this additivity is not guaranteed in the second model. Theoretical arguments were 
put forward both for and against long-range energy transfer, 21 ’ and attempts were made 
to carry out the forcing experiment, but they were hindered by the low Reynolds numbers 
available from simulations at the time. The matter has remained controversial until now, 
and our additive data on Fig. 3 favors the local model. 

We now focus on the scaling of At. We start by introducing the strong assumption that 
the values of At we measure are between scales r within an inertial range, where r and e are 
the only relevant quantities. Within this simplistic framework, then, a velocity in r-space 
can be defined as 


r = e/p B , (4) 

where the energy density p (r) is a real-space equivalent of E{k). We put forward the 
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FIG. 4. (a) Derivative of q = UiUi/2 with respect to r. The dashed line corresponds to e 2 ^ 3 r 1 ^ 3 . 
(b) Energy content within a band of scales between r and r + A r, where A r goes from a given r 
to the next bigger r in the plotted series. A q = 4 UiUi(r + Ar) — 4 UiUi(r ). 


following candidate for p E , based on q = mui/2: 


Pe 


dq 

dr 


( 5 ) 


Other densities have been introduced in physical space. Townsend’s r-derivative of the 
correlation function, or the signature function found in Ref. 24 are two examples. We 
chose our expression as it is based on the filtering approach we use. The data on Fig. 4(a) 
shows that the energy density —dq/dr in our two flows is a positive quantity. In Fig. 4(b) 
we see that the energy —A q contained between r and r + A?’ is a quantity which grows 
proportionaly to r 2 / 3 within a reasonable range in HIT3 - not so in HST with much smaller 
scale separation. Such r 2 / 3 behaviour is consistent with the Kolmogorov-Obukhov theory, 2,8 
yet it is based on a completely different flow decomposition from the structure function. 


FIG. 5. E 1/3 At r ^r a /r 2 a 13 
and (E(r a )), with r > r a . 



Symbols as in Table III. The solid line corresponds to Eq. (6) 


Substituting p E in Eq. (4) by e 2 ^ 3 r 1//3 , and integrating from r a to r leads to 

A^ ra =iT 1 / 3 (r 2 / 3 -r2/ 3 ), (6) 

which implies that our data should fall on a straight line when plotted logarithmically 
as done in Fig. 5. The agreement is not completely unsatisfactory. Particularly when 
considering the inertial range assumptions used which have no reason to apply if r or r a are 
below the viscous limit of 30?7, or given the poor compliance of HST to the inertial range 
scaling with r - see Fig. 4. A dashed line following Eq. (6) for HIT3 was added on Fig. 2(c). 

A derivation of Eq. (6) carried out in spectral space can be found in Ref. 12, leading to a 
k~ 2 / 3 dependence of At. Ref. 25 studied Lagrangian time correlations of S(fc) and £, which 
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hinted at a k~ 2 / 3 dependence of the peaks in their correlations - see inset of their Fig. 2. 
Earlier, the same group measured the temporal correlation between the energy at a given 
scale and the energy at a smaller scale found later by following the flow in both forced and 
decaying HIT. Their conclusion that the peak in correlation happens later for increasing 
scale separation is consistent with what we observe. A difficulty with that work was the use 
of correlations of energy rather than energy flux. Fluxes are the quantities conserved across 
cascades, and the natural objects for their study. Furthermore, they could only consider 
one-jump delays, ruling out the additivity test on Fig. 3 which supports the locality of the 
energy cascade in an average sense. 
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